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Abstract 

Saturation affects a significant rate of images recorded by the Atmospheric 
Imaging Assembly on the Solar Dynamics Observatory. This paper describes 
a computational method and a technological pipeline for the de-saturation 
of such images, based on several mathematical ingredients like Expectation 
Maximization, image correlation and interpolation. An analysis of the com¬ 
putational properties and demands of the pipeline, together with an assess¬ 
ment of its reliability are performed against a set of data recorded from the 
Feburary 25 2014 flaring event. 
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1. Introduction 

As in any optical system, telescopes introduce blur into the images they 
produce. The characteristics of such a blur are encoded in the instrument 
point spread function (PSF), which describes the response of the whole opti¬ 
cal hardware to the signal produced by a distant point source. Deconvolving 
any blurring and artifacts associated with this PSF is a crucial step among 
the variety of tasks to be addressed in image processing. Unlike many instru¬ 
ments, whose PSF is made of just a main central peak associated with signal 
diffusion, some telescopes have a PSF with a double structure, whereby a sec¬ 
ond, more complex component overlays the diffusion core, containing peaks 
that replicate the central one according to regular patterns of varying inten¬ 
sity EU- Such replications are caused by diffraction which is produced by 
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the waves scattering off the regular wire meshes that support the instrument 
entrance filters. 

During intense solar flares the diffraction fringes appear as artefacts at 
long range from a central core but these are a timely source of information 
on the true image as the direct information has been destroyed by saturation 
effects in the recording hardware. In fact, the principle of CCD devices [6J is 
to convert photons into electrons. The electrons that accumulate in the CCD 
wells are read out and converted to data numbers (DNs) by software. The 
amount of charge that a single pixel can accumulate is finite and coarsely 
related to its area. However, the probability of trapping an electron within a 
pixel well decreases when the well is approaching saturation, which implies 
that the linear relation between light intensity and signal degrades and the 
response of a saturated pixel drops. For intense photon fluxes, the major 
fraction of the signal accumulates in the CCD image pixels only up to the 
saturation limit while another part is coherently scattered by the support 
to produce the diffraction fringes in the remote pixels without any deviation 
from a linear response. Therefore, the diffraction signal is unaffected by 
saturation; further, the diffraction pattern is a direct signature of the photon 
flux in the saturated region and therefore, in principle, an inverse process 
can restore the original scene. 

Image saturation has been an issue for a number of telescopes in solar as¬ 
tronomy, like the NASA Solar Terrestrial Relations Observatory (STEREO) 
[9] and the NASA Transition Region and Coronal Explorer (TRACE) [4J. 
However, with the introduction of the Atmospheric Imaging Assembly (AIA) 
in the Solar Dynamics Observatory (SDO) [5j de-saturation has become a 
big data issue. Indeed AIA provides full-disk 4096 x 4096 pixel images of the 
solar chromosphere and corona in seven EUV wavelengths with a time ca¬ 
dence of 12 seconds and saturation affects a large fraction of these images: a 
coarse estimate, based on the fraction of saturated frames in typical C-class, 
M-class and X-class events and on the typical number of such events provided 
by the RHESSI flare-list, suggests that more than 10 5 AIA images per year 
are saturated. This pathology prevents a systematic exploitation of AIA in 
both the accomplishment of its scientific objectives and the integration with 
data from other observation missions. 

DESAT is a software package now available in Solar Software (SSW) 
that realizes an automatic de-saturation of AIA images by using a correla¬ 
tion/inversion analysis of the diffraction fringes produced in the telescope ob¬ 
servations [7]. The conceptual features of DESAT are essentially two. First, 
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the tool strongly relics on the shape and characteristics of the instrument 
PSF in both its components, the diffusion core and the diffraction PSF. Sec¬ 
ond, the mathematical content at the basis of DESAT is not at all trivial and 
involves Expectation Maximization (EM) [8j for realizing image segmenta¬ 
tion and ” inverse” diffraction. Moreover, this same mathematical framework 
makes DESAT easily applicable to image de-saturation in the case of instru¬ 
ments other then SDO/AIA , under the condition that these instruments are 
characterized by a similar two-component PSF. 

DESAT pipeline also provides a module for reducing secondary effects 
typical of instruments like AIA. When a pixel reaches saturation, additional 
charge cannot be accommodated in the pixel well and therefore spills over 
into adjacent pixels, often causing them to saturate. This spread of charge 
to a pixel region around the primary saturation core is named blooming and 
typically appears as bright ray along a preferential direction depending on the 
CCD design as well as saturated pixels surrounding the primary saturation. 
DESAT takes care of blooming and ameliorates its effects by exploiting the 
fact that each saturated AIA image is embedded in a time sequence of variable 
duration, made of acquisitions performed with smaller exposure times and 
therefore non-saturated. DESAT utilizes the information contained in the 
sequence of non-saturated images to interpolate an estimate of the pixel 
content in the blooming region. Although the applications considered so far 
show that this approach is rather reliable, we note that it is significantly less 
robust than the segment of DESAT pipeline reducing primary saturation: 
the routines of such segment, in fact, utilize information provided by the 
instrument at the right time and directly related to primary saturation by 
the physics of diffraction. 

The plan of the paper is as follows. Section 2 provides the mathematical 
setup for DESAT. Section 3 describes the pipeline and Section 4 performs 
some numerical tests to validate the effectiveness of the toolbox. Our con¬ 
clusions are offered in Section 5, together with a discussion of the main open 
issues. 

2. Mathematical setup 

This Section briefly summarizes the mathematical theory at the basis 
of the de-saturation approach described in much more details in [lOj. The 
starting point for this theory is that the inverse diffraction problem at the 
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Figure 1: An example of saturated AIA image with highlighted the overall (primary + 
blooming) saturation region and the diffraction fringes. The event is the one occurred on 
February 25 2014, time point 00:45:56 UT. 


basis of the DESAT toolbox can be described by the linear model 

g = Kf + b, (1) 

where (see Figure [I]) g is the recorded signal in correspondence of the diffrac¬ 
tion fringes, / is the true flux in the primary saturation region, K is the 
diffraction component of the AIA PSF mapping the primary saturation re¬ 
gion onto the diffraction pattern and b is the portion of the image background 
associated to the diffraction fringes. 

Then the mathematical challenge of solving equation (JT]) is many-fold and 
requires the following steps: 

Segmentation. This step separates the image to de-saturate in four re¬ 
gions: the one containing the diffraction fringes, the primary saturation 
region, the blooming region and the remaining part of the image, which 
will stay untouched during the overall de-saturation process. We first of 
all observe that the diffraction pattern and the primary saturation re¬ 
gion are directly connected, since the diffraction fringes are significantly 
illuminated, via the diffraction PSF, just by the flux directed toward 
the pixels contained in the primary saturation region. The crucial is¬ 
sue of this step is therefore to segment the saturation and the blooming 
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regions. In fact, at a first inspection of the image, such regions cannot 
be distinguished and are made by those pixels whose grey level con¬ 
tent is bigger than 16,000 DN pixel -1 , i.e. the saturation threshold 
for AIA. However, the diffraction PSF strongly correlates just with the 
pixel in the primary saturation, while the correlation between K and 
the blooming part of the region must be very weak. Therefore DESAT 
computes the correlation map by means of one iteration of the EM al¬ 
gorithm and retains as primary saturation pixels just those pixels for 
which such map (after projection onto the data space by means of the 
diffusion PSF) has components larger than the saturation value. 

Reconstruction. This step solves equation (JTJ) by means of EM. Assuming 
that the diffraction fringes g are affected by Poisson noise (CCD data 
are more correctly quasi-Poisson but the approximation is reasonable 
in this context), EM implements maximum likelihood under positivity 
constraint by optimally stopping the iteration [H] 

fk+l T f _ 9 _\ 

1 K T 1 \Kf k + b) ’ 1 j 

where 1 is made of all ones. Several different stopping rules can be 
applied for EM. We empirically found that the one mathematically 
discussed in |2] and applied in [3] for image reconstruction from count 
modulation profiles recorded by the Reuven Ramaty High Energy Solar 
Spectroscopic Imager (RHESSI) works very appropriately also in the 
case of AIA data and therefore we used it in DESAT. 

Image synthesis. The map provided by EM in the reconstruction step 
is an estimate of the incoming flux in correspondence to the primary 
saturation region. Since saturation occurs in the data region, DESAT 
projects the EM solution onto the data space by means of the core 
diffusion PSF and then synthesize the de-saturated image by glueing 
the result of the projection onto the image background. 

All three steps described above exploit the knowledge of an estimate of 
the image background. DESAT constructs this background map via interpo¬ 
lation of the pixel content in non-saturated images acquired in a time range 
containing the time point when the saturated image has been recorded. Once 
the background map is available different sub-regions of it are used by DE¬ 
SAT for different purposes: the sub-region corresponding to the diffraction 
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fringes is used in the iterative solution of equation ([!]) by means of EM (it 
is matrix b in (§ and g); the subregion where the overall saturation region 
(primary saturation + blooming) is projected by diffraction is used to com¬ 
pute correlation in the segmentation process; the subregion corresponding to 
the blooming region is used in the synthesis of the final de-saturated image. 

3. The DESAT pipeline 

The pipeline realized by DESAT in Interactive Data Language (IDL) and 
currently available in SSW is illustrated in Figure [2] and is made of the 
following steps: 

INPUT: this step loads a hie containing a set of saturated and unsaturated 
AIA images, corresponding to a selected time interval and a vector 
indicating the specific wavelengths to process. These two ingredients 
are saved in the ’save’ step. 

read: this step simply reads the hies in the datasets associated to the 
selected wavelengths. 

background estimation: this is one of the computational cores of the 
pipeline. It provides an estimate of the background for the whole held 
of view covered by the image. 

PSF computation: this step computes the diffusion PSF component, the 
diffraction PSF component and the global PSF. 

segmentation: this step utilizes correlation to identify the position of 
primary saturation pixels, of blooming pixels and of the pixels in the 
diffraction fringes. 

EM reconstruction: this step utilizes EM to reconstruct the incoming 
EUV hux in correspondence of the primary saturation region. 

image synthesis: this step projects the reconstructed hux in the image 
space by means of the diffusion PSF and glues the result on the back¬ 
ground to provide the unsaturated image. 

OUTPUT: the output of the pipeline is made of a dataset of de-saturated 
AIA images which are saved again in the ’save’ step. 
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The most mathematical aspects of the pipeline are described in detail in 
the boxes ’Algorithm 1’ and ’Algorithm 2’, respectively. 

Specifically, step 1 in the ’Algorithm 1’ box de-convolves all non-saturated 
images in order to reduce the effects of diffraction and diffusion. Then, steps 
from 2 through 6 construct a corresponding set of filtered Fourier trans¬ 
formed images (a Butterworth low-pass filter is also applied and, for each 
transformed image, just pixels where the filter values are greater than a 
threshold are retained). For each pixel, steps 7 through 11 construct two 
time sequences (for the real and imaginary parts of the Fourier transformed 
images) such that each time point contains the pixel intensity of the corre¬ 
sponding image in the sequence. Steps 12 and 14 use quadratic interpolation 
to infer the pixel values at those time points corresponding to the saturated 
images. Steps from 15 through 21 re-arrange these interpolation values to 
construct Fourier transform maps of the background at the saturation time 
points (in particular, steps 19 through 21 fill up the image regions cast away 
by the Butterworth filter by using the pixel contents of the Fourier transform 
of a reference unsaturated image). Steps 19 through 21 perform the inverse 
Fourier transform of the background images and finally step 24 projects such 
images onto the data space by means of the core diffusion PSF. We point 
out that working in the Fourier domain notably decreases the computational 
burden of the whole process, since most information in these data is con¬ 
centrated in a limited region of small frequencies, where the interpolation is 
performed. 

The ’Algorithm 2’ box describes the pipeline routines that identify the 
primary saturation region with respect to the blooming one, apply EM to 
restore the flux in the primary saturation region and project this information 
into the data space in order to synthesize the de-saturated image. Specifically, 
steps from 2 through 9 realize the segmentation of the saturated region by 
identifying the primary saturation region and determine the location of the 
diffraction fringes. This is obtained by thresholding the correlation between 
the diffraction PSF and the set of pixels contaminated by diffraction under 
the assumption that the whole saturation region were primary (the correla¬ 
tion map is computed by one iteration of EM). Step 10 through 12 restore 
the photon flux in the primary saturation region by numerically solving equa¬ 
tion (|TJ) via EM. Eventually, steps 13 through 25 provide the de-saturated 
image synthesis by glueing together the de-saturated information and the 
background. 
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Figure 2: The DESAT pipeline. 
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4. Numerical tests 

We tested DESAT in the case of the flaring event of February 25 2014, 
in the time interval between 00:45:00 UT and 00:47:00 UT, which presented 
nine saturated images with different saturation degree (i.e., different numbers 
of saturated pixels). An example of such images, corresponding to the time 
point 00:45:56, is represented in Figure[3j which is the same one as in Figure]!] 
the right panel in the same figure shows the same image after the application 
of DESAT. In order to assess the relative and overall computational cost 
of the pipeline and its experimental reliability, we generated the saturated 
images according to three different sizes and tracked the computational cost 
of the ’segmentation’, ’EM reconstruction’ and ’image synthesis’ elements of 
the pipeline. The results of this test are in Table [l] for the images with size 
500 x 500, Table [2] for the images with size 750 x 750 and Table [3] for the 
images with size 1000 x 1000, respectively. These results show that the biggest 
impact on the overall computational time of the pipeline comes from the EM 
step; this is expected, since this step requires several matrix times product 
computation. On the other hand the segmentation and image synthesis steps 
rather mildly depend on both the number of saturated/data pixel and the 
dimension of the image under de-saturation. 

As far as the reliability of the process is concerned, we computed the C- 
statistic values that measure the discrepancy between the observed data in 
the region of the diffraction fringes and the expectation values predicted by 
the reconstructed flux. The last column of the three Tables reports the aver¬ 
aged value of these C-statistic values and show that it significantly increases 
with the number of saturation points to cure. Further, Table [4] shows the be¬ 
havior of the averaged C-statistic values when computed in correspondence of 
the fringe at positions with increasing distance from the image center (we re¬ 
ported here just the results concerning the 500 x 500 size, but the other cases 
behave coherently). Interestingly, the reliability of the de-saturation process 
seems not to be significantly influenced by the position of the reproduced 
diffraction fringes with respect to the de-saturated core. 

5. Conclusions 

DESAT is a first step to the full exploitation of SDO/AIA images: the 
availability of a fully automatic technological pipeline will make EUV data 
available also in the case of highly energetic flaring events (although sat¬ 
uration sometimes occurs even in the case of mildly energetic flares), and 
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Figure 3: DESAT at work. Left panel: the same saturated image as in Figure [T] Right 
panel: the result of de-saturation. 


time 

# sat 

# data 

seg 

EM (time) 

EM (iter) 

synthesis 

C-stat 

00:45:08 

1304 

112426 

12.1 

46.6 

58 

10.7 

34.8 

00:45:23 

563 

95712 

12.7 

14.6 

53 

10.8 

2.3 

00:45:32 

1459 

110570 

11.5 

64.0 

72 

10.6 

46.8 

00:45:48 

23 

16337 

11.8 

2.2 

37 

10.4 

3.9 

00:45:56 

1592 

121232 

11.7 

75.4 

88 

10.6 

57.9 

00:46:12 

26 

20097 

11.7 

3.0 

42 

10.3 

5.9 

00:46:20 

1981 

45691 

11.8 

64.4 

71 

10.1 

26.4 

00:46: 36 

147 

45691 

11.8 

6.4 

43 

10.2 

5.5 

00:46:44 

2300 

108488 

11.4 

71.4 

67 

10.2 

21.6 


Table 1: DESAT pipeline at work against images with 500 x 500 size. First column: 
recording time. Second column: number of saturated pixels. Third column: number of 
pixels in the diffraction fringe. Fourth column: computational time for the ’segmentation’ 
step. Fifth column: computational time for the ’EM reconstruction’ step. Sixth column: 
number of iteration for EM. Seventh column: computational time for the ’image synthesis’ 
step. Eighth column: averaged C-statistic values computed in correspondence with the 
diffraction fringes. 
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time 

# sat 

# data 

seg 

EM (time) 

EM (iter) 

synthesis 

C-stat 

00:45:08 

1304 

195097 

15.5 

280.1 

55 

10.7 

39.0 

00:45:23 

563 

153760 

15.8 

25.8 

50 

108 

2.0 

00:45:32 

1459 

193362 

15.8 

599.8 

67 

10.2 

29.4 

00:45:48 

23 

20911 

17.5 

2.2 

42 

10.4 

2.0 

00:45:56 

1592 

201723 

18.7 

502.0 

92 

10.3 

36.1 

00:46:12 

26 

26327 

15.9 

20.0 

37 

10.4 

20.0 

00:46:20 

1981 

205386 

16.2 

502.0 

72 

10.3 

26.9 

00:46:36 

147 

62513 

15.9 

2.0 

42 

10.0 

2.4 

00:46:44 

2300 

210492 

16.3 

368.1 

73 

10.4 

13.4 


Table 2: DESAT pipeline at work against images with 750 x 750 size. First column: 
recording time. Second column: number of saturated pixels. Third column: number of 
pixels in the diffraction fringe. Fourth column: computational time for the ’segmentation’ 
step. Fifth column: computational time for the ’EM reconstruction’ step. Sixth column: 
number of iteration for EM. Seventh column: computational time for the ’image synthesis’ 
step. Eighth column: averaged C-statistic values computed in correspondence with the 
diffraction fringes. 


time 

# sat 

# data 

seg 

EM (time) 

EM (iter) 

synthesis 

C-stat 

00:45:08 

1304 

237279 

18.3 

282.5 

55 

8.3 

10.9 

00:45:23 

563 

178711 

17.0 

46.1 

49 

8.1 

1.6 

00:45:32 

1459 

239203 

16.3 

356.0 

69 

8.4 

36.1 

00:45:48 

23 

21312 

16.3 

2.2 

26 

8.2 

1.6 

00:45:56 

1592 

246624 

16.6 

485.4 

94 

8.5 

21.1 

00:46:12 

26 

28784 

16.9 

2.2 

36 

8.3 

20.0 

00:46:20 

1981 

254602 

16.6 

401.8 

75 

8.5 

21.8 

00:46:36 

147 

72714 

18.3 

5.5 

38 

8.8 

1.7 

00:46:44 

2300 

263798 

21.5 

404.2 

78 

8.6 

12.4 


Table 3: DESAT pipeline at work against images with 1000 x 1000 size. First column: 
recording time. Second column: number of saturated pixels. Third column: number of 
pixels in the diffraction fringe. Fourth column: computational time for the ’segmentation’ 
step. Fifth column: computational time for the ’EM reconstruction’ step. Sixth column: 
number of iteration for EM. Seventh column: computational time for the ’image synthesis’ 
step. Eighth column: averaged C-statistic values computed in correspondence with the 
diffraction fringes. 


11 





25-Feb-2014 

Ri 

R-2 

Rs 

00:45:08 

2.21037 

2.47712 

2.37171 

00:45:23 

0.991982 

1.39406 

1.16723 

00:45:32 

2.36914 

2.76948 

2.54095 

00:45:48 

1.49663 

2.11597 

1.74499 

00:45:56 

7.87822 

10.2115 

9.05733 

00:46:12 

1.93773 

1.68562 

1.81033 

00:46:20 

6.91328 

4.70813 

5.73614 

00:46:36 

2.37601 

1.69629 

1.82247 

00:46:44 

4.32539 

4.28034 

4.01359 


Table 4: Averaged C-statistic values computed in correspondence of the diffraction fringes 
in three image portions at increasing distance from the saturated core. R\ corresponds to 
the region closest to the core (between 21 arcsec and 60 arcsec from the image center). i ?2 
corresponds to the intermediate region (between 60 arcsec to 105 arcsec). R 3 corresponds 
to the most peripheral region (between 105 and 150 arcsec). 


will realize this information improvement in the framework of a ’big data’ 
approach, which is coherent with the SDO original concept. 

However, DESAT has still limitations and therefore it will probably pave 
the way to a series of actions that will ameliorate its performance at different 
levels. For example, there exists a reliability issue that should be addressed 
in DESAT and which is testified by still high C-statistic values, mainly in 
the case of strongly saturated data. To reduce such values, three aspects 
should be considered. First, all methods in DESAT strongly rely on the 
form of the PSF applied and therefore one should study the impact of the 
model of the PSF (particularly, its diffraction component) on the quality 
of de-saturation (DESAT is currently using the synthetic estimate of the 
diffraction PSF provided by SSW). Second, AIA filters provide data that are 
centered around a specific wavelength but that contain information coming 
from a range of wavelengths around the central one: accounting for this 
wavelength-dependent dispersion effect (which involves also the instrument 
PSF) should most likely help to reduce the C-statistic values. Third, the way 
DESAT is currently estimating the image background and therefore the way 
information is restored in the blooming region is most likely non-optimal and 
other strategies should be explored. 

Finally, at a more computational level, a significant improvement should 
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come from the implementation of an automatic recipe to select the size of 
the image to de-saturate, in such a way that the corresponding input data 
will represent the optimal trade-off between the information content in the 
fringes and the computational burden. 
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Algorithm 1 background 

Require: A dataset of AIA saturated and not saturated images corre¬ 
sponding to a fixed time interval. Let {s(p,j), p = 1, • • •, N pixe i, 
j = 0, • • •, N sat — 1} be the set of saturated images and {ns(p, i), 
p = 1, • • •, Npixei i = 0, • • •, N no _ sat — 1} the set of unsaturated images. 
Let U be the set of pixels in each image. Let PSF and PSF C be the 
complete PSF and the core PSF of the AIA instrument respectively. 
Ensure: A dataset of images {bg(-,j), j = 0, • • • ,N sat — 1} providing back¬ 
ground estimation for saturated images based on information from not 
saturated ones. 

1: Deconvolve images ns(-, •) by using PSF: ns '(•, •) y- ns(-, •) 

2: for all i = 0,..., N no sat - 1 do 

3: Compute the Fourier Transform of image ns'(-, i ): Fns'(-, i ) y- ns'(i ) 

4: Apply a fourth-order Butterworth low-pass filter to image Fns'(-,i ): 

BFns'(-,i) y- Fns'(-,i ) 

5: end for 

6: Identify the set S C U of pixels where the filter values are greater than 
a threshold value 
7: for all p G S do 
8: for alii = 0,, N no _ sat - 1 do 

9: Fill entry i of a vector reaLorig: reaForig(i ) y- iR.(BFns'(p,i )) 

10: Fill entry i of a vector imag-orig: imagjorig{i) y- ^ (BFns'(p,i)) 

11: end for 

12: Apply interpolation to vector reaLorig: realJnterp y- reaLorig 

13: Apply interpolation to vector imag-orig: imagJnterp y- imag.orig 

14: end for 

15: for all j = 0,..., N sat - 1 do 

16: for all p E S do 

17: Nbg'{j>,j) y- comp\ex(real-interp(j), imagJnterp(j)) 

18: end for 

19: for all p G U \ S do 

20: Fbg'(p,j) y- Fns'ipS)) 

21: end for 

22: Compute the Inverse Fourier Transform of image Fbg'(-,j): bg'(-,j) y- 

Nbg'(;j) 

23: end for 

24: Convolve images bg'(-,-) by using PSFq: bg(-,-) y- bg'(- } •) 
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Algorithm 2 segmentation - EM reconstruction - image synthesis 
Require: A dataset of AIA saturated images {s(p,j), p = 1, • • •, N pixe i, 
j = 0, • • •, N sat — 1} and corresponding background estimations {bg(p, j), 
p = 1, • • •, Npi xe i , j = 0, • • •, N sat — 1}. Let U be the set of pixels in each 
image. Let PSF, PSF C and PSF D be the complete PSF, the core PSF 
and the diffraction PSF of the AIA instrument respectively. 

Ensure: A dataset of desaturated images { s desat (j), J = 0,--- ,N sat -l}. 
1: for all j = 0, ..., N sat - 1 do 
2 : Identify the set S' C U of saturated pixels 

3: Identify the set MS — U \ S' of unsaturated pixels 

4: Convolve the indicator function I5/ with PSFd to identify the set 

r C MS of pixels where diffraction from pixels in S' virtually occurs: 

r <- i s > * psf d 

5: Compute the correlation C between the image values in F' and PSFp 

given the background estimation in T' 

6: Convolve C with PSFc to obtain a correlation map C in the image 

space: C C' * PSF C 

7: Threshold C to identify the set S C S' of primary saturation pixels 

8: Identify the set B C S' of bloomed pixels: B <— S' \S 

9: Convolve I5 with PSFd to identify the set T C T’ of pixels where 

diffraction from pixels in S actually occurs: T I 5 * PSFo 
10 : Restrict s(-,j) image to F: sj?(j) s(p,j), p G F 

11 : Restrict bg(-,j ) image to F\ bgj?(j ) bg(p, j), p e F 

12 : Apply EM algorithm to solve the inverse problem given by Eq. (fTI) : 

xAj) <r- EM(sjr(j), bgjr(j), PSF D ) 

13: Restrict bg(-,j ) image to B: bgjs(j) bg(p, j), p G B 

14: for all p G U do 

15: if p G F then 

16: S desat (p,j ) «- [x F (j) * PSF c ](p ) 

17: else if p G B then 

18: s desat (p,j) «- bg B (p,j) 

19: else if p G F then 

20: S desat (p,j ) <- s(p,j ) - [x F (j) * PSF D }(p ) 

21: else 

22: S desat (p,j) 4- s(p,j) 

23: end if 

24: end for 

25: end for 
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